All code and data for this workshop is available at the following URL:

https://github.com/kaybenleroll/data_workshops

Code is available in the ws_pgm_202009/ directory.

Content in this workshop is based on the book Probabalistic Graphical Models: Principles and Techniques by Soren Hojsgaard.

Also look at the vignettes for the packages gRain and gRbase

Remember that this topic is massive. I could easily give a full semester course on this stuff to really do it justice, so most of this workshop is just me working through the material as I learn it.

As a result, it is highly likely this worksheet and code contains typos, errors, logical flaws and other mistakes in need of correction in this workshop, so if you note any, please let me know so I can try to fix them!

If you want to look into this topic more, there is an old Coursera course by Daphne Koller (tough going but excellent):

https://www.coursera.org/learn/probabilistic-graphical-models

This course was based on her textbook

Probabalistic Graphical Models: Principles and Techniques

The Gaussian Graphical Models section of the workshop is based on material found in the book “Bayes Networks With Examples in R” by Scutari and Denis.

1 Basic Concepts

A graph is a mathematical object that can be defined as a pair

\[ \mathcal{G} = (V, E), \]

where \(V\) is a set of vertices or nodes, and \(E\) is a set of edges that joins two vertices. Edges in general may be directed, undirected or bidirected. They are typically visualised by using shapes or points for the nodes and lines for the edges.

The concept of conditional independence is related to that of statistical independence. Suppose we have three random variables \(A\), \(B\) and \(C\), then \(A\) and \(B\) are conditionally independent given \(C\), written

\[ A \perp B \mid C, \]

iff, for every given value \(c\) in \(C\), \(A\) and \(B\) are independent in the conditional distribution given \(C = c\).

Another way of saying this is that for some \(f\) a generic density or probability mass function, then one characteristic of \(A \perp B \, | \, C\) is that

\[ f(a, b \mid c) = f(a \mid c) f(b \mid c). \]

An equivalent characterisation is that the joint density of \(A\), \(B\) and \(C\) factorises as

\[ f(a, b, c) = g(a, c) \, h(b, c). \]

Finally, we will also make heavy use of Bayes’ Rule, the standard formula for relating conditional probabilities:

\[ P(A \mid B) = \frac{P(A, B)}{P(B)} = \frac{P(B \mid A) P(A)}{P(B)}. \]

1.1 Conditional Probability

To explain the concept of conditional probability we first look at some basic exercises in probability.

In the following questions, we are rolling two six-sided dice (denoted 2D6) with \(D_1\) and \(D_2\) denoting the result of the first and second die respectively, and we denote \(T\) as being the sum of the two - i.e.

\[ T = D_1 + D_2 \]

Question: We start simple - with no further detail, what is the probability of the total being 11?

Answer: We do not know anything about \(D_1\) or \(D_2\), so there are two outcomes that result in \(T = 11\):

\[(D_1 = 5, D_2 = 6) \text{ or } (D_1 = 6, D_2 = 5).\]

Thus, the probability of rolling 11 is

\[ P(T = 11) = \frac{2}{36} = 0.05556 \]


Question: What is the probability of getting 11 if the first dice is 5?

Answer: Out total is 11 if the second dice is 6, so there is only 1 possible outcome of those six possible rolls.

\[ P(T = 11 \mid D_1 = 5) = \frac{1}{6} = 0.1667 \]

1.1.1 Exercises

We now ask a few questions to ensure you understand these concepts.

  1. If the second dice is a 4, what is the probability the total was 7?
  2. What is the probability of the total being 10 or greater?
  3. If we know the first dice roll is 4, what is the probability of the total being 10 or greater?
  4. If we know the second dice roll is 2, what is the probability of the total being 10 or greater?
  5. What is the probability of the total being 3 or less?

1.2 Conditional Dependence

Now suppose we know the totals of each dice instead? Suppose we know that \(T = 9\) and we have the following distribution for \(D_1\):

\[ P(D_1) = \{0, \, 0, \, 0, \, 0, \, 0.5, \, 0.5 \} \] induces the following probability distribution for \(D_2\):

\[ P(D_2) = \{0, \, 0, \, 0.5, \, 0.5, \, 0, \, 0 \} \]

because \(D_1 = 5 \implies D_2 = 4\) and \(D_1 = 6 \implies D_2 = 3\)

Thus, while we conceptually think of \(D_1\) and \(D_2\) as being independent of one another, they become dependent if we have knowledge of \(T\).

In this case we say that \(D_1\) and \(D_2\) are conditionally dependent on one another given \(T\).

Now we define new variables, \(X_1\) and \(X_2\):

\[ X_1 = \begin{cases} 1 \text{ iff } T \text{ is even}\\ 0 \text{ otherwise} \end{cases} \; \; \; X_2 = \begin{cases} 1 \text{ iff } T >= 9\\ 0 \text{ otherwise} \end{cases} \]

Given no further information, we say that \(X_1\) and \(X_2\) are dependent, as knowing information about one variables affects our knowledge of the other.

However, what happens if we know something about \(T\)?

In this case, we say that \(X_1\) and \(X_2\) are independent given \(T\).

2 Bayesian Network

We start with Bayesian Networks - where the nodes on the graph represent discrete random variables.

2.1 The Sprinkler Network

We start with a basic example of a Bayes Network: the sprinkler network.

yn <- c("Yes", "No")

cptable_lst <- list(
  cptable(~Rain,                        levels = yn, values = c(2, 8)),
  cptable(~Sprinkler + Rain,            levels = yn, values = c(1, 99, 4, 6)),
  cptable(~wetGrass + Rain + Sprinkler, levels = yn, values = c(99, 1, 8, 2, 9, 1, 0, 1))
)

sprinkler_cptlist <- cptable_lst %>% compileCPT()
sprinkler_grain   <- sprinkler_cptlist %>% grain()

sprinkler_grain %>% plot()

Two events can cause grass to be wet: Either the sprinkler is on or it is raining. Rain has a direct effect on the use of the sprinkler (namely that when it rains, the sprinkler is usually not turned on).

This can be modeled with a Bayesian network. The variables (R)ain, (S)prinkler, Wet(G)rass have two possible values: (y)es and (n)o.

We can factorise the joint probability mass function as

\[ p_{GSR}(g, s, r) = p_{G \mid SR}(g \mid s, r) p_{S \mid R}(s \mid r) p_R(r) \]

or overloading the notation a little:

\[ P(G, S, R) = P(G \mid S, R) \; P(S, R) = P(G \mid S, R) \; P(S \mid R) \; P(R) \]

This means we can construct the joint probability table by starting with the \(conditional probability tables\) (CPTs).

Create the 3 CPTs using the parray() function and the following conditional probabilities:

\[\begin{align*} P(R) &= 0.2 \\ P(S|R) &= 0.01 & P(S|\neg R) &= 0.4\\ P(G|S,R) &= 0.99 & P(G|S, \neg R) &= 0.9 & P(G|\neg S, R) &= 0.8 & P(G|\neg S, \neg R) &= 0 \end{align*}\]

First we want to construct this network using our various conditional probabilities:

yn <- c("yes", "no")

## P(R)
p_R <- parray(
  varNames = "Rain",
  levels = list(yn),
  values = c(0.2, 0.8)
)

## P(S|R)
p_S_R <- parray(
  varNames = c("Sprinkler", "Rain"),
  levels   = list(yn, yn), 
  values   = c(0.01, 0.99, 0.4, 0.6)
)

## P(G|S,R)
p_G_SR <- parray(
  varNames = c("GrassWet", "Sprinkler", "Rain"),
  levels   = list(yn, yn, yn),
  values   = c(0.99, 0.01, 0.8, 0.2, 0.9, 0.1, 0, 1)
)

ftable(p_G_SR, row.vars = "GrassWet")
##          Sprinkler  yes        no     
##          Rain       yes   no  yes   no
## GrassWet                              
## yes                0.99 0.90 0.80 0.00
## no                 0.01 0.10 0.20 1.00

We now combine these probabilities into a full joint distribution p_GSR:

p_GSR <- tabListMult(
  list(p_G_SR, p_S_R, p_R)
)

ftable(p_GSR, row.vars = "GrassWet")
##          Sprinkler     yes              no        
##          Rain          yes      no     yes      no
## GrassWet                                          
## yes                0.00198 0.28800 0.15840 0.00000
## no                 0.00002 0.03200 0.03960 0.48000

No suppose we know that the grass is wet. What is the probability it is raining?

p_RG  <- tabMarg(p_GSR, c("Rain", "GrassWet"))  ## P(R,G)
p_G   <- tabMarg(p_RG, "GrassWet")              ## P(G)
p_R_G <- tabDiv(p_RG, p_G)                      ## P(R|G)

Reading across this CPT, we see that if the grass is wet, there is about a 36% chance it is due to rain.

While the above methods of manipulating CPTs works, the package gRain provides us better functionality for answering questions such as this given observed evidence.

3 Genetic Inheritance

We now turn our attention to analysing genetic inheritance on the chromosomes for a given DNA sequence.

An allele is the DNA sequence at a marker and can take two values marked \(A\) or \(B\) (in practice there can be 10 or 20 different values).

A genotype is an unordered pair of alleles: \(AA\), \(AB\), or \(BB\).

The genotype of a person at a specific marker is a random variable with state space \(\{AA, AB, BB\}\).

We are interested in the joint distribution of genotypes for a group of people.

We need to make a number of assumptions regarding genetic inheritance that we list here:

  • A child inherits one allele from each parent independently.
  • The parent’s two alleles have equal probability of being passed on to the child.
  • Each combination has probability \(0.25\); some lead to the same genotype for the child.

So in this case we have the the joint probability distribution as being

\[ P(m, f, c) = P(m) \, P(f) \, P(c \mid m, f) \]

## # A tibble: 27 x 4
##    child mother father  prob
##    <chr> <chr>  <chr>  <dbl>
##  1 AA    AA     AA      1   
##  2 AA    AA     AB      0.5 
##  3 AA    AA     BB      0   
##  4 AA    AB     AA      0.5 
##  5 AA    AB     AB      0.25
##  6 AA    AB     BB      0   
##  7 AA    BB     AA      0   
##  8 AA    BB     AB      0   
##  9 AA    BB     BB      0   
## 10 AB    AA     AA      0   
## # … with 17 more rows
## Rows: 27
## Columns: 4
## $ child  <chr> "AA", "AA", "AA", "AA", "AA", "AA", "AA", "AA", "AA", "AB", "A…
## $ mother <chr> "AA", "AA", "AA", "AB", "AB", "AB", "BB", "BB", "BB", "AA", "A…
## $ father <chr> "AA", "AB", "BB", "AA", "AB", "BB", "AA", "AB", "BB", "AA", "A…
## $ prob   <dbl> 1.00, 0.50, 0.00, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…

Suppose we have a population frequency of alleles being 70% \(A\) and 30% \(B\). We want to calculate the distribution of genotypes in the population.

This is a straightforward application of the Binomial distribution, and we use the R function dbinom()

genotype_probs <- dbinom(0:2, size = 2, prob = 0.3)

print(genotype_probs)
## [1] 0.49 0.42 0.09

Thus we have it that 49% of the population have genotype \(AA\), 42% are \(AB\) and 9% are \(BB\).

We now want to construct the CPTs for each of the nodes on the network:

mother_cpt <- cptable(~ mother, values = genotype_probs, levels = genotypes)
print(mother_cpt)
## {v, pa(v)} :
##  chr "mother"
##    [,1]
## AA 0.49
## AB 0.42
## BB 0.09
father_cpt <- cptable(~ father, values = genotype_probs, levels = genotypes)
print(father_cpt)
## {v, pa(v)} :
##  chr "father"
##    [,1]
## AA 0.49
## AB 0.42
## BB 0.09
p_inheritance <- prob_tbl %>%
  arrange(father, mother, child) %>%
  pull(prob)

child_cpt <- cptable(
  ~ child | mother + father,
  values = p_inheritance,
  levels = genotypes
  )

print(child_cpt)
## {v, pa(v)} :
##  chr [1:3] "child" "mother" "father"
##    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## AA    1  0.5    0  0.5 0.25  0.0    0  0.0    0
## AB    0  0.5    1  0.5 0.50  0.5    1  0.5    0
## BB    0  0.0    0  0.0 0.25  0.5    0  0.5    1

Using the CPTs, we construct the Bayesian network:

genetic_family_grain <- list(child_cpt, mother_cpt, father_cpt) %>%
  compileCPT() %>%
  grain()

plot(genetic_family_grain)

3.1 Querying the Network

By using the gRain package it allows us to use the in-built functionality to query the network.

To start with, we want to see the marginal distribution of alleles for the father.

genetic_family_grain %>%
  querygrain(nodes = "father")
## $father
## father
##   AA   AB   BB 
## 0.49 0.42 0.09

We can also look at the joint distribution of the mother and child:

genetic_family_grain %>%
  querygrain(
    nodes = c("child", "mother"),
    type  = "joint"
  )
##      mother
## child    AA    AB    BB
##    AA 0.343 0.147 0.000
##    AB 0.147 0.210 0.063
##    BB 0.000 0.063 0.027

We can pass this output to ftable() if we want:

genetic_family_grain %>%
  querygrain(
    nodes = c("child", "mother"),
    type  = "joint"
    ) %>%
  ftable(col.vars = "child")
##        child    AA    AB    BB
## mother                        
## AA           0.343 0.147 0.000
## AB           0.147 0.210 0.063
## BB           0.000 0.063 0.027

We can also produce the conditional distributions for each of the possible values for the father given values of the mother and child.

genetic_family_grain %>%
  querygrain(
    nodes = c("father", "child", "mother"),
    type  = "conditional"
    ) %>%
  ftable(col.vars = "father")
##              father   AA   AB   BB
## child mother                      
## AA    AA            0.70 0.30 0.00
##       AB            0.70 0.30 0.00
##       BB             NaN  NaN  NaN
## AB    AA            0.00 0.70 0.30
##       AB            0.49 0.42 0.09
##       BB            0.70 0.30 0.00
## BB    AA             NaN  NaN  NaN
##       AB            0.00 0.70 0.30
##       BB            0.00 0.70 0.30

3.2 Paternity Test

Now suppose we know that a child has genotype \(AB\), and the mother has genotype \(BB\). Given that a man has genotype \(AB\), what can we say about the chances of the man being the father of the child?

This is a more complicated question than it sounds, as we need to think about a few different ideas. We make use of the pEvidence() and setEvidence().

p_fmc <- genetic_family_grain %>%
  setEvidence(
    evidence = list(mother = "BB", child = "AB", father = "AB")
    )

p_fmc %>% print()
## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "child" "mother" "father"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 mother             TRUE         BB
## 2  child             TRUE         AB
## 3 father             TRUE         AB
##   pEvidence: 0.018900
p_fmc %>% pEvidence()
## [1] 0.0189

So we now have a probability for finding this particular combination of alleles in a mother / father / child grouping.

So now we calculate the probability of any given man being \(AB\)?

p_f <- genetic_family_grain %>%
  setEvidence(evidence = list(father = "AB"))

p_f %>% print()
## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "child" "mother" "father"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 father             TRUE         AB
##   pEvidence: 0.420000
p_f %>% pEvidence()
## [1] 0.42

Finally, we look at the probability of the child being \(AB\) with a mother \(BB\).

p_mc <- genetic_family_grain %>%
  setEvidence(
    evidence = list(mother = "BB", child = "AB")
    )

p_mc %>% print()
## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:3] "child" "mother" "father"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 mother             TRUE         BB
## 2  child             TRUE         AB
##   pEvidence: 0.063000
p_mc %>% pEvidence()
## [1] 0.063

So now we need to look at the various ratios of these probabilities to assess probabilities.

p_man <- pEvidence(p_fmc) / pEvidence(p_f)
print(p_man)
## [1] 0.045
p_father <- p_man / pEvidence(p_mc)
print(p_father)
## [1] 0.7142857

3.3 Extended Genetic Family

We can now extend this model to add extended families, such as uncle/aunts and grandparents.

c_mf <- parray(
  varNames = c("child", "mother", "father"),
  levels   = rep(list(genotypes), 3),
  values   = p_inheritance
  )

f_gmgf <- parray(
  varNames = c("father", "grandmother", "grandfather"),
  levels   = rep(list(genotypes), 3),
  values = p_inheritance
  )

u_gmgf <- parray(
  varNames = c("uncle", "grandmother", "grandfather"),
  levels   = rep(list(genotypes), 3),
  values   = p_inheritance
  )

m  <- parray("mother",      levels = list(genotypes), values = genotype_probs)
gf <- parray("grandfather", levels = list(genotypes), values = genotype_probs)
gm <- parray("grandmother", levels = list(genotypes), values = genotype_probs)


extended_family_grain <- list(c_mf, m, f_gmgf, u_gmgf, gm, gf) %>%
  compileCPT() %>%
  grain()

extended_family_grain %>% iplot()

Now suppose the man is dead and we do not have any genetic information for him, but we know his brother tests as \(AA\). How does this affect the evidence for potential paternity?

p_uncle <- extended_family_grain %>%
  setEvidence(evidence = list(mother = "BB", child = "AB", uncle = "AA"))

p_uncle %>% print()
## Independence network: Compiled: TRUE Propagated: TRUE 
##   Nodes: chr [1:6] "child" "mother" "father" "uncle" "grandmother" "grandfather"
##   Evidence:
##    nodes is.hard.evidence hard.state
## 1 mother             TRUE         BB
## 2  child             TRUE         AB
## 3  uncle             TRUE         AA
##   pEvidence: 0.037485
p_uncle %>% pEvidence()
## [1] 0.037485

We also want to show how we calculate similar values by using querygrain().

extended_family_grain %>%
  querygrain(
    nodes = c("child", "mother", "uncle"),
    type  = "joint"
    ) %>%
  ftable(col.vars = "uncle")
##              uncle       AA       AB       BB
## child mother                                 
## AA    AA           0.204085 0.123480 0.015435
##       AB           0.087465 0.052920 0.006615
##       BB           0.000000 0.000000 0.000000
## AB    AA           0.036015 0.082320 0.028665
##       AB           0.102900 0.088200 0.018900
##       BB           0.037485 0.022680 0.002835
## BB    AA           0.000000 0.000000 0.000000
##       AB           0.015435 0.035280 0.012285
##       BB           0.006615 0.015120 0.005265

4 Chest-Clinic Network

We now turn our attention to the “Chest-Clinic Network” as an example of a Bayesian network.

Based on research by Lauritzen and Spiegelhalter, we construct a network based on the following domain knowledge:

Shortness-of-breath (dyspnoea) may be due to tuberculosis, lung cancer or bronchitis, or none of them, or more than one of them. A recent visit to Asia increases the chances of tuberculosis, while smoking is known to be a risk factor for both lung cancer and bronchitis. The results of a single chest X-ray do not discriminate between lung cancer and tuberculosis, as neither does the presence or absence of dyspnoea.

We use this knowledge to construct a Bayesian network - breaking up the above into discrete pieces of knowledge.

chestclinic_dag <- list(
    "asia",
    c("tub", "asia"),
    "smoke",
    c("lung", "smoke"),
    c("bronc", "smoke"),
    c("either", "lung", "tub"),
    c("xray", "either"),
    c("dysp", "bronc", "either")
    ) %>%
  dag()

chestclinic_dag %>% plot()

4.1 Constructing the Bayesian Network

a    <- cptable(~ asia,
                values = c(1, 99),
                levels = yn)
t_a  <- cptable(~ tub | asia,
                values = c(5, 95, 1, 99),
                levels = yn)
s    <- cptable(~ smoke,
                values = c(5, 5),
                levels = yn)
l_s  <- cptable(~ lung | smoke,
                values = c(1, 9, 1, 99),
                levels = yn)
b_s  <- cptable(~ bronc | smoke,
                values = c(6, 4, 3, 7),
                levels = yn)
e_lt <- cptable(~ either | lung + tub,
                values = c(1, 0, 1, 0, 1, 0, 0, 1),
                levels = yn)
x_e  <- cptable(~ xray | either,
                values = c(98, 2, 5, 95),
                levels = yn)
d_be <- cptable(~ dysp | bronc + either,
                values = c(9, 1, 7, 3, 8, 2, 1, 9),
                levels = yn)

chestclinic_cpt_grain <- list(a, t_a, s, l_s, b_s, e_lt, x_e, d_be) %>%
  compileCPT() %>%
  grain()

chestclinic_cpt_grain %>% summary()
## Independence network: Compiled: TRUE Propagated: FALSE 
##  Nodes : chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" "xray" "dysp"
##  Number of cliques:                 6 
##  Maximal clique size:               3 
##  Maximal state space in cliques:    8

We first should check some of the various conditional probabilities to ensure we have built the network correctly:

chestclinic_cpt_grain %>%
  querygrain(
    nodes = c("xray", "either"),
    type  = "conditional"
  )
##      either
## xray   yes   no
##   yes 0.98 0.05
##   no  0.02 0.95
chestclinic_cpt_grain %>%
  querygrain(
    nodes = c("dysp", "bronc", "either"),
    type  = "conditional"
    ) %>%
  ftable(row.vars = c("either", "bronc"))
##              dysp yes  no
## either bronc             
## yes    yes        0.9 0.1
##        no         0.7 0.3
## no     yes        0.8 0.2
##        no         0.1 0.9

4.2 Building Networks from Data

Now that we have our network, we can use data to construct the CPTs. Assuming no data is missing, we can use a tibble of data to set the probabilities.

data(chestSim500)
data(chestSim1000)
data(chestSim10000)


chestdata_500_tbl   <- chestSim500   %>% as_tibble()
chestdata_1000_tbl  <- chestSim1000  %>% as_tibble()
chestdata_10000_tbl <- chestSim10000 %>% as_tibble()

chestdata_500_tbl   %>% glimpse()
## Rows: 500
## Columns: 8
## $ asia   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ tub    <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ smoke  <fct> no, no, yes, yes, yes, no, yes, no, no, no, no, no, yes, no, y…
## $ lung   <fct> no, no, no, no, no, no, no, no, no, no, no, no, yes, no, no, n…
## $ bronc  <fct> yes, yes, yes, yes, yes, no, yes, no, no, yes, no, no, yes, ye…
## $ either <fct> no, no, no, no, no, no, no, no, no, no, no, no, yes, no, no, n…
## $ xray   <fct> no, no, no, yes, no, no, no, no, no, no, no, no, yes, no, no, …
## $ dysp   <fct> no, yes, yes, no, yes, yes, yes, no, no, yes, no, no, yes, yes…
chestdata_1000_tbl  %>% glimpse()
## Rows: 1,000
## Columns: 8
## $ asia   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ tub    <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ smoke  <fct> no, yes, yes, no, yes, yes, yes, no, yes, no, no, no, yes, yes…
## $ lung   <fct> no, no, no, no, no, yes, no, no, no, no, no, no, no, no, no, n…
## $ bronc  <fct> yes, yes, no, no, yes, yes, yes, no, no, no, no, no, yes, yes,…
## $ either <fct> no, no, no, no, no, yes, no, no, no, no, no, no, no, no, no, n…
## $ xray   <fct> no, no, no, no, no, yes, no, no, no, no, no, no, no, no, no, n…
## $ dysp   <fct> yes, yes, no, no, yes, yes, yes, no, no, no, no, no, yes, yes,…
chestdata_10000_tbl %>% glimpse()
## Rows: 10,000
## Columns: 8
## $ asia   <fct> no, no, no, no, no, no, no, no, no, yes, yes, no, no, no, no, …
## $ tub    <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ smoke  <fct> no, yes, yes, yes, no, no, yes, yes, yes, no, no, yes, yes, ye…
## $ lung   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ bronc  <fct> no, yes, no, yes, no, yes, no, yes, yes, yes, yes, no, no, no,…
## $ either <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ xray   <fct> no, no, no, no, no, no, no, no, no, no, no, no, no, no, no, no…
## $ dysp   <fct> no, yes, no, yes, no, yes, no, yes, no, yes, yes, no, no, no, …

We now can build multiple networks for each of the datasets we have.

chestsim_500_grain <- chestclinic_dag %>%
  grain(data = chestdata_500_tbl, smooth = 0.1) %>%
  compile(propagate = TRUE)

chestsim_1000_grain <- chestclinic_dag %>%
  grain(data = chestdata_1000_tbl, smooth = 0.1) %>%
  compile(propagate = TRUE)

chestsim_10000_grain <- chestclinic_dag %>%
  grain(data = chestdata_10000_tbl, smooth = 0.1) %>%
  compile(propagate = TRUE)

The unconditional probability of a patient having lung cancer should match the proportion of cases in the dataset.

chestsim_500_grain %>%
  querygrain(nodes = "lung", type = "joint")
## lung
##        yes         no 
## 0.04636229 0.95363771
chestdata_500_tbl %>%
  count(lung, name = "count") %>%
  mutate(prob = count / sum(count))
## # A tibble: 2 x 3
##   lung  count  prob
##   <fct> <int> <dbl>
## 1 yes      23 0.046
## 2 no      477 0.954

We do the same thing for the other networks we constructed.

First we look at the 1,000 row dataset.

chestsim_1000_grain %>%
  querygrain(nodes = "lung", type = "joint")
## lung
##        yes         no 
## 0.05317813 0.94682187
chestdata_1000_tbl %>%
  count(lung, name = "count") %>%
  mutate(prob = count / sum(count))
## # A tibble: 2 x 3
##   lung  count  prob
##   <fct> <int> <dbl>
## 1 yes      53 0.053
## 2 no      947 0.947

Now we look at the network based on 10,000 datapoints.

chestsim_10000_grain %>%
  querygrain(nodes = "lung", type = "joint")
## lung
##        yes         no 
## 0.05571778 0.94428222
chestdata_10000_tbl %>%
  count(lung, name = "count") %>%
  mutate(prob = count / sum(count))
## # A tibble: 2 x 3
##   lung  count   prob
##   <fct> <int>  <dbl>
## 1 yes     557 0.0557
## 2 no     9443 0.944

4.3 Using Evidence

We can use given evidence to adjust the conditional probabilities.

Suppose the individual has visited Asia and has dyspnoea, what is the conditional probability that the person has lung cancer?

chestsim_500_grain %>%
  querygrain(
    nodes    = "lung",
    evidence = list(asia = "yes", dysp = "yes"),
    type     = "marginal"
    )
## $lung
## lung
##        yes         no 
## 0.06414233 0.93585767

Given this evidence we see the probability of lung cancer has now increased to around 6.5%.

We repeat this exercise for the 1,000 network and see what probability it shows.

chestsim_1000_grain %>%
  querygrain(
    nodes    = "lung",
    evidence = list(asia = "yes", dysp = "yes"),
    type     = "marginal"
    )
## $lung
## lung
##       yes        no 
## 0.1030665 0.8969335

In this network, the probability is 10%.

chestsim_10000_grain %>%
  querygrain(
    nodes    = "lung",
    evidence = list(asia = "yes", dysp = "yes"),
    type     = "marginal"
    )
## $lung
## lung
##      yes       no 
## 0.103872 0.896128

As before, the probability is about 10%.

4.4 Marginal Proportions for Diseases

We also want to calculate the different marginal probabilities for the three conditions - Tuberculosis, Lung cancer and Dyspnoea - according to the different

chestclinic_cpt_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )
## $tub
## tub
##    yes     no 
## 0.0104 0.9896 
## 
## $lung
## lung
##   yes    no 
## 0.055 0.945 
## 
## $bronc
## bronc
##  yes   no 
## 0.45 0.55
chestsim_500_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )
## $tub
## tub
##        yes         no 
## 0.01432307 0.98567693 
## 
## $lung
## lung
##        yes         no 
## 0.04636229 0.95363771 
## 
## $bronc
## bronc
##       yes        no 
## 0.4540341 0.5459659
chestsim_1000_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )
## $tub
## tub
##         yes          no 
## 0.007197136 0.992802864 
## 
## $lung
## lung
##        yes         no 
## 0.05317813 0.94682187 
## 
## $bronc
## bronc
##       yes        no 
## 0.4360235 0.5639765
chestsim_10000_grain %>%
  querygrain(
    nodes = c("lung", "bronc", "tub"),
    type  = "marginal"
    )
## $tub
## tub
##        yes         no 
## 0.01091936 0.98908064 
## 
## $lung
## lung
##        yes         no 
## 0.05571778 0.94428222 
## 
## $bronc
## bronc
##       yes        no 
## 0.4422023 0.5577977

5 Additional Topics

Before we finish the topic of Bayesian networks with categorical variables, there are a number of topics worth discussion briefly: in particular, how to scale the model to larger numbers of variables and levels, and how to use the data to determine possible configurations of the network in terms of relationships between variables.

5.1 Scaling Up Bayesian Networks

All of the above approaches are example of the `Brute Force’ approach which is done by calculating the full joint distribution for the network \(p(V)\) as a multiple of the CPTs that comprise it.

\[ p(V) = p(a) \, p(t \mid a) \, p(s) \, p(l \mid s) \, p(b \mid s) \, p(e \mid t, l) \, p(d \mid e, b) \, p(x \mid e) \]

This gives \(p(V)\) represented by a table with \(2^8 = 256\) entries.

We then marginalise and condition as desired to calculate whatever probabilities we need.

The scaling of this approach is bad. A network with 80 variables, each with 10 values has a joint probability space of \(10^{80}\), approximately the count of atoms in the universe.

We are going to need a bigger boat…

So, we must find an approach that does not require the full joint distribution, instead focusing on the the low dimensional CPTs and send `messages’ between them.

To use a network it first needs to be compiled and then propagated.

Compilation of a network based on CPTs is first moralized — edges are added between the parents of each node, and then directed edges are replaced with undirected ones. It is then triangulated to form a triangulated graph.

The CPTs are transformed into clique potentials defined on the cliques of the chordal graph.

We see this process as below:

chestclinic_moralized    <- chestclinic_dag %>% moralize()
chestclinic_triangulated <- chestclinic_moralized %>% triangulate()

chestclinic_triangulated %>% plot()

Once we create the DAG, we view the triangulated data as a junction tree. The messages passed between connected nodes on the graph involve the common variables in the nodes, and propagating the information involves a double pass down and up the tree.

We now repeat the process from the above section, but now use the this approach to calculate our probabilities.

chestsim500_ev_grain <- chestsim_500_grain %>%
  setFinding(nodes = "asia", states = "yes", propagate = FALSE) %>%
  setFinding(nodes = "dysp", states = "yes", propagate = FALSE) %>%
  propagate()

We now use this network to perform the same calculations. As a reminder, this is calculating the marginal, joint and conditional probability for lung cancer and bronchitis given that the person recently visited Asia and exhibited shortness of breath (dyspnoea).

chestsim500_ev_grain %>%
  querygrain(
    nodes = c("lung", "bronc"),
    type  = "marginal"
    )
## $lung
## lung
##        yes         no 
## 0.06414233 0.93585767 
## 
## $bronc
## bronc
##       yes        no 
## 0.6761692 0.3238308
chestsim500_ev_grain %>%
  querygrain(
    nodes = c("lung", "bronc"),
    type  = "joint"
    )
##      bronc
## lung         yes         no
##   yes 0.03832565 0.02581668
##   no  0.63784353 0.29801414
chestsim500_ev_grain %>%
  querygrain(
    nodes = c("lung", "bronc"),
    type  = "conditional"
    )
##      bronc
## lung         yes         no
##   yes 0.05668056 0.07972275
##   no  0.94331944 0.92027725

5.2 Learning Network Structure

One final topic to think about is using the data to infer a model to explain it - so called model selection or structural learning.

Note that rather than using the data to determine the parameters of the model, we instead use the data to infer the structure of the model, that is, what the Bayesian network will look like.

For the chest clinic, we already know what the structure was, so we can compare the output of the various structural learning approaches to see how well it does.

Before we do that, we will plot the ‘true’ network:

chestclinic_dag %>% plot()

We now try a number of approaches to learn the structure - starting with the hill climbing algorithm implemented via hc().

chestSim500 %>%
  hc() %>%
  amat() %>%
  ggm::essentialGraph() %>%
  as("graphNEL") %>%
  plot()

We see that most of the network structure has been ‘learned’ from the data, apart from the relationship between “asia” and “tub” - presumably as this effect is small and thus weakly signalled within the data.

We repeat this exercise with chestSim10000 as this may allow us to pick that relationship out.

chestSim10000 %>%
  hc() %>%
  amat() %>%
  ggm::essentialGraph() %>%
  as("graphNEL") %>%
  plot()

We can try one larger set of data with 100,000 entries.

data(chestSim100000)

chestSim100000 %>%
  hc() %>%
  amat() %>%
  ggm::essentialGraph() %>%
  as("graphNEL") %>%
  plot()

6 Networks with Continuous Variables

All our Bayesian networks up to this point have focused on networks with discrete variables, and we want to extend our knowledge to building and using these networks.

For the purposes of this workshop, we will focus on Gaussian Bayesian networks: networks where all the continuous variables are distributed according to Gaussian densities (both univariate and multivariate) as well as the following assumptions:

  • every node follows a normal distribution,
  • nodes without any parent, known as root nodes, are described by the respective marginal distributions,
  • the conditioning effect of the parent nodes is given by an additive linear term in the mean, and does not affect the variance,
  • the local distribution of each node can be equivalently expressed as a Gaussian linear model which includes an intercept and the node’s parents, without any interaction terms.

6.1 Crop Analysis

We start with an example where we wish to analyse the crop yield of a particular plant. To do this, we have a number of things to consider:

  1. The genetic and environmental potential of the plant
  2. Vegatative mass
  3. Harvest mass - the crop

For the first point, we construct two latent variables, \(G\) and \(E\) to measure the (G)enetic potential and (E)nvironmental potential of the plant.

We also want to measure various aspects of the (V)egatative mass, as these are indicators to the quality of the plant and hence affect the quality and quantity of the crop yield. Thus, we also have the variable \(V\).

The crop is a function of the vegatative mass and so is measured by the (N)umber of seeds and their mean (W)eight - giving us variables \(N\) and \(W\).

Finally, the (C)rop is determined by the number of seeds and the weight.

Converting all of the above into a Bayesian network, we can visualise it as follows:

Plot of the Crop Network

(#fig:construct_crop_bn_dag)Plot of the Crop Network

If we were to specify a full joint probability distribution for \(p\) variables, this means we need \(p\) means, \(p\) variances and \(\frac{1}{2}p(p-1)\) elements of the correlation matrix and need to ensure that the correlation matrix is positive definite.

Thus, this problem scales as \(\mathcal{O}(N^2)\).

By using this graphical model, we only need to specify the various local, conditional distributions, simplifying the problem.

6.2 Constructing the Network

Before we add these distributions we want to construct the model. We can specify the above conditional factorisation of the model with the function model2network() function from the bnlearn package.

plant_crop_dagbnlearn <- model2network("[G][E][V|G:E][N|V][W|V][C|N:W]")

plant_crop_dagbnlearn %>% print()
## 
##   Random/Generated Bayesian network
## 
##   model:
##    [E][G][V|E:G][N|V][W|V][C|N:W] 
##   nodes:                                 6 
##   arcs:                                  6 
##     undirected arcs:                     0 
##     directed arcs:                       6 
##   average markov blanket size:           2.67 
##   average neighbourhood size:            2.00 
##   average branching factor:              1.00 
## 
##   generation algorithm:                  Empty

Suppose for example, we specify the following conditional distributions:

\[\begin{eqnarray*} G &\sim& \mathcal{N}(50, 10^2) \\ E &\sim& \mathcal{N}(50, 10^2) \\ V \mid G = g, E = e &\sim& \mathcal{N}(-10.35534 + 0.5 g + 0.70711 e, 5^2) \\ N \mid V = v &\sim& \mathcal{N}(45 + 0.1 v, 9.949874 ^ 2) \\ W \mid V = v &\sim& \mathcal{N}(15 + 0.7 v, 7.141428 ^ 2) \\ C \mid N = n, W = w &\sim& \mathcal{N}(0.3 n + 0.7 w, 6.25 ^ 2) \end{eqnarray*}\]

We now set up these conditional distributions in bnlearn. From inspection we see we are setting the parameters of a regression model at each node. As discussed, each node is distributed according to a Gaussian density, with conditional dependence manifesting as the conditional mean being a linear function of the dependent variables.

dist_E <- list(coef = c("(Intercept)" = 50), sd = 10)
dist_G <- list(coef = c("(Intercept)" = 50), sd = 10)
dist_V <- list(coef = c("(Intercept)" = -10.35534, E = 0.70711, G = 0.5), sd = 5)
dist_N <- list(coef = c("(Intercept)" = 45, V = 0.1), sd = 9.949874)
dist_W <- list(coef = c("(Intercept)" = 15, V = 0.7), sd = 7.141428)
dist_C <- list(coef = c("(Intercept)" = 0, N = 0.3, W = 0.7), sd = 6.25)

distrib_lst = list(E = dist_E, G = dist_G, V = dist_V, N = dist_N, W = dist_W, C = dist_C)


plant_crop_bnlearn <- custom.fit(plant_crop_dagbnlearn, dist = distrib_lst)

6.3 Inference on the Model

We now have constructed this ‘toy’ GBN and wish to focus on how we use this network and model to answer questions about the system - much like we did before in the discrete case.

Due to the use of Gaussian densities, we can use exact inference on this network. This has the benefit of calculating exact conditional probabilities, but imposes limits on the types of questions that can be ask.

As a result, we will focus on methods of approximate inference - that is, the use of Monte Carlo simulation to provide statistical estimates of the parameters of interest.

This gives us much more flexibility in terms of the types of questions we can ask, at the cost of producing answers that are only approximately correct.

That said, the error in this estimate is more often than not much smaller than any errors due to the choice of the model.

6.3.1 Naive Simulation

We start with a simple, but naive method of calculating these conditional expectations: we create a sample of \(N\) iterations from the realisations of the joint distribution, and use this sample to calculation estimates of the expectations.

n_sim <- 1000000

plant_crop_sample_tbl <- rbn(n = n_sim, plant_crop_bnlearn) %>% as_tibble()

plant_crop_sample_tbl %>% glimpse()
## Rows: 1,000,000
## Columns: 6
## $ C <dbl> 44.80154, 57.24168, 65.23836, 62.43541, 53.70927, 64.48227, 50.8389…
## $ E <dbl> 57.81968, 41.31241, 48.04928, 55.16027, 53.30161, 58.30459, 64.3550…
## $ G <dbl> 49.37286, 63.04870, 72.86645, 36.11139, 47.21211, 48.66679, 56.3595…
## $ N <dbl> 48.00326, 71.30687, 60.71959, 48.65638, 46.91299, 45.58467, 57.6481…
## $ V <dbl> 52.10191, 49.71834, 66.85698, 55.27992, 56.06012, 56.21025, 58.1256…
## $ W <dbl> 55.09897, 57.57701, 62.78116, 67.73014, 56.41255, 60.74187, 42.8630…

6.3.1.1 Marginal Distribution of C

We start with some simple questions: what is the marginal distribution of \(C\)?

ggplot(plant_crop_sample_tbl) +
  geom_histogram(aes(x = C), binwidth = 1) +
  scale_y_continuous(labels = label_comma()) +
  xlab("Value of C") +
  ylab("Simulation Count") +
  ggtitle("Marginal Distribution of C")

This is exactly what we would expect - we see the marginal distribution of \(C\) follows that of a Gaussian density.

6.3.1.2 E[C | E = 50]

Now suppose we now that the crop was grown in an environment with average potential - that is \(E = 50\). What is the conditional distribution of \(C\)?

To do this, we filter the distribution for \(E = 50\) and look at the resulting values.

plant_crop_sample_tbl %>% filter(E == 50)
## # A tibble: 0 x 6
## # … with 6 variables: C <dbl>, E <dbl>, G <dbl>, N <dbl>, V <dbl>, W <dbl>

This approach needs to be amended, as we see that none of the samples has \(E = 50\) exactly. This is not an accident, for continuous distributions, the probability mass of any single value is 0. Instead, we need to look at a narrow range of values close to the value we care about.

plant_crop_e50_5_tbl    <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 5.0)
plant_crop_e50_1_tbl    <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 1.0)
plant_crop_e50_01_tbl   <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 0.1)
plant_crop_e50_001_tbl  <- plant_crop_sample_tbl %>% filter(abs(E - 50) < 0.01)

plant_crop_e50_tbl <- list(
    `Width 5.0`  = plant_crop_e50_5_tbl,
    `Width 1.0`  = plant_crop_e50_1_tbl,
    `Width 0.1`  = plant_crop_e50_01_tbl,
    `Width 0.01` = plant_crop_e50_001_tbl
    ) %>%
  bind_rows(.id = "width_label")

plot_tbl <- plant_crop_e50_tbl %>%
  group_by(width_label) %>%
  summarise(
    .groups = "drop",
    
    mean_value = mean(C),
    
    p10 = quantile(C, 0.10),
    p25 = quantile(C, 0.25),
    p50 = quantile(C, 0.50),
    p75 = quantile(C, 0.75),
    p90 = quantile(C, 0.90)
    )

ggplot(plot_tbl) +
  geom_errorbar(aes(x = width_label, ymin = p25, ymax = p75),
                width = 0, size = 3, colour = "red") +
  geom_errorbar(aes(x = width_label, ymin = p10, ymax = p90),
                width = 0, size = 1, colour = "black") +
  geom_point(aes(x = width_label, y = p50),        colour = "red") +
  geom_point(aes(x = width_label, y = mean_value), colour = "black") +
  expand_limits(y = 0) +
  xlab("Interval Width") +
  ylab("Value of C") +
  ggtitle("Conditional Values of C")

We can create a more direct comparison using facetted histograms:

ggplot(plant_crop_e50_tbl) +
  geom_histogram(aes(x = C), binwidth = 1) +
  facet_wrap(vars(width_label), scales = "free_y", ncol = 2) +
  xlab("C") +
  ylab("Frequency Count") +
  ggtitle("Facetted Histogram Plot")

This visualisation could be better so we instead try to use a density plot.

ggplot(plant_crop_e50_tbl) +
  geom_line(aes(x = C, colour = width_label), stat = "density") +
  xlab("C") +
  ylab("Density") +
  ggtitle("Comparison Plot of Estimated Density Functions")

6.3.1.3 E[C | E = 70]

Now we repeat all of the above for \(E = 70\).

plant_crop_e70_tbl <- list(
    `Width 5.0`  = plant_crop_sample_tbl %>% filter(abs(E - 70) < 5.0),
    `Width 1.0`  = plant_crop_sample_tbl %>% filter(abs(E - 70) < 1.0),
    `Width 0.1`  = plant_crop_sample_tbl %>% filter(abs(E - 70) < 0.1),
    `Width 0.01` = plant_crop_sample_tbl %>% filter(abs(E - 70) < 0.01)
    ) %>%
  bind_rows(.id = "width_label")

plant_crop_e70_tbl %>% glimpse()
## Rows: 72,093
## Columns: 7
## $ width_label <chr> "Width 5.0", "Width 5.0", "Width 5.0", "Width 5.0", "Widt…
## $ C           <dbl> 55.16532, 65.78010, 53.84612, 52.69393, 56.45567, 61.6004…
## $ E           <dbl> 68.56438, 67.05995, 71.12048, 66.88804, 71.49582, 69.3356…
## $ G           <dbl> 54.32818, 56.42899, 52.84883, 55.80996, 41.14224, 41.3920…
## $ N           <dbl> 52.35925, 57.40292, 50.67462, 60.10386, 59.95746, 43.6093…
## $ V           <dbl> 67.22984, 67.47149, 68.79282, 57.69178, 62.78158, 54.2236…
## $ W           <dbl> 61.06274, 56.41881, 66.60216, 59.49617, 52.83071, 63.2764…
ggplot(plant_crop_e70_tbl) +
  geom_line(aes(x = C, colour = width_label), stat = "density") +
  xlab("C") +
  ylab("Density") +
  ggtitle("Comparison Plot of Estimated Density Functions")

6.3.1.4 P(E, G | C > 80)

Now suppose we want to look at the environmental and genetic values associated with very good crops (say \(C > 80\)).

plant_crop_goodcrop_tbl <- plant_crop_sample_tbl %>%
  filter(C > 80)

ggplot(plant_crop_goodcrop_tbl) +
  geom_point(aes(x = E, y = G)) +
  xlab("E") +
  ylab("G") +
  ggtitle("Conditional Joint Distribution of E and G")

6.3.2 Importance Sampling

While the simulation method works for the more commonly-occuring scenarios, there is an issue for rare combinations of events as most of the generated iterations will not be relevant for the analysis and the computation work is wasted.

We use the cpdist() and cpquery() functions to allow us to do this work without fully realising each of the iterations.

plant_crop_importance_tbl <- cpdist(
    plant_crop_bnlearn,
    nodes    = c("E", "G"),
    evidence = (C > 80),
    n = 1000000
    ) %>%
  as_tibble()

plant_crop_importance_tbl %>% glimpse()
## Rows: 1,342
## Columns: 2
## $ E <dbl> 47.46442, 68.75327, 67.12545, 60.86390, 48.37759, 59.48258, 67.9352…
## $ G <dbl> 60.10531, 60.65128, 48.64428, 53.58544, 55.32707, 67.23726, 78.7154…

Now that we have generated the samples, we can plot both \(E\) and \(G\) variables.

ggplot(plant_crop_importance_tbl) +
  geom_point(aes(x = E, y = G)) +
  xlab("Environmental Potential") +
  ylab("Genetic Potential") +
  ggtitle("Plot of High Crop Yield Environmental Variables")

We can now compare these two scatterplots.

comparison_tbl <- list(
    simulation = plant_crop_goodcrop_tbl,
    importance = plant_crop_importance_tbl
    ) %>%
  bind_rows(.id = "source")

ggplot(comparison_tbl) +
  geom_point(aes(x = E, y = G, colour = source)) +
  xlab("Environmental Potential") +
  ylab("Genetic Potential") +
  ggtitle("Comparison Plot of (E,G) Where C > 80")

For further use of Bayesian Networks, the following packages are likely of interest:

7 R Environment

sessioninfo::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.2 (2020-06-22)
##  os       Ubuntu 20.04.1 LTS          
##  system   x86_64, linux-gnu           
##  ui       RStudio                     
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Etc/UTC                     
##  date     2020-12-01                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  assertthat     0.2.1    2019-03-21 [1] RSPM (R 4.0.0)
##  backports      1.1.10   2020-09-15 [1] RSPM (R 4.0.2)
##  BiocGenerics   0.36.0   2020-10-27 [1] Bioconductor  
##  blob           1.2.1    2020-01-20 [1] RSPM (R 4.0.0)
##  bnlearn      * 4.6.1    2020-09-21 [1] RSPM (R 4.0.2)
##  bookdown       0.20     2020-06-23 [1] RSPM (R 4.0.2)
##  broom          0.7.1    2020-10-02 [1] RSPM (R 4.0.2)
##  cellranger     1.1.0    2016-07-27 [1] RSPM (R 4.0.0)
##  cli            2.1.0    2020-10-12 [1] RSPM (R 4.0.2)
##  colorspace     1.4-1    2019-03-18 [1] RSPM (R 4.0.0)
##  conflicted   * 1.0.4    2019-06-21 [1] RSPM (R 4.0.0)
##  cowplot      * 1.1.0    2020-09-08 [1] RSPM (R 4.0.2)
##  crayon         1.3.4    2017-09-16 [1] RSPM (R 4.0.0)
##  DBI            1.1.0    2019-12-15 [1] RSPM (R 4.0.0)
##  dbplyr         1.4.4    2020-05-27 [1] RSPM (R 4.0.0)
##  digest         0.6.25   2020-02-23 [1] RSPM (R 4.0.0)
##  dplyr        * 1.0.2    2020-08-18 [1] RSPM (R 4.0.2)
##  ellipsis       0.3.1    2020-05-15 [1] RSPM (R 4.0.0)
##  evaluate       0.14     2019-05-28 [1] RSPM (R 4.0.0)
##  fansi          0.4.1    2020-01-08 [1] RSPM (R 4.0.0)
##  farver         2.0.3    2020-01-16 [1] RSPM (R 4.0.0)
##  forcats      * 0.5.0    2020-03-01 [1] RSPM (R 4.0.0)
##  fs             1.5.0    2020-07-31 [1] RSPM (R 4.0.2)
##  generics       0.0.2    2018-11-29 [1] RSPM (R 4.0.0)
##  ggm            2.5      2020-02-16 [1] RSPM (R 4.0.2)
##  ggplot2      * 3.3.2    2020-06-19 [1] RSPM (R 4.0.1)
##  glue           1.4.2    2020-08-27 [1] RSPM (R 4.0.2)
##  gRain        * 1.3-6    2020-07-30 [1] RSPM (R 4.0.2)
##  graph          1.68.0   2020-10-27 [1] Bioconductor  
##  gRbase       * 1.8-6.7  2020-07-03 [1] RSPM (R 4.0.2)
##  gtable         0.3.0    2019-03-25 [1] RSPM (R 4.0.0)
##  haven          2.3.1    2020-06-01 [1] RSPM (R 4.0.2)
##  highr          0.8      2019-03-20 [1] RSPM (R 4.0.0)
##  hms            0.5.3    2020-01-08 [1] RSPM (R 4.0.0)
##  htmltools      0.5.0    2020-06-16 [1] RSPM (R 4.0.1)
##  httr           1.4.2    2020-07-20 [1] RSPM (R 4.0.2)
##  igraph         1.2.6    2020-10-06 [1] RSPM (R 4.0.2)
##  jsonlite       1.7.1    2020-09-07 [1] RSPM (R 4.0.2)
##  knitr          1.30     2020-09-22 [1] RSPM (R 4.0.2)
##  labeling       0.3      2014-08-23 [1] RSPM (R 4.0.0)
##  lattice        0.20-41  2020-04-02 [2] CRAN (R 4.0.2)
##  lifecycle      0.2.0    2020-03-06 [1] RSPM (R 4.0.0)
##  lubridate      1.7.9    2020-06-08 [1] RSPM (R 4.0.2)
##  magrittr     * 1.5      2014-11-22 [1] RSPM (R 4.0.0)
##  MASS           7.3-51.6 2020-04-26 [2] CRAN (R 4.0.2)
##  Matrix         1.2-18   2019-11-27 [2] CRAN (R 4.0.2)
##  memoise        1.1.0    2017-04-21 [1] RSPM (R 4.0.0)
##  modelr         0.1.8    2020-05-19 [1] RSPM (R 4.0.0)
##  munsell        0.5.0    2018-06-12 [1] RSPM (R 4.0.0)
##  pillar         1.4.6    2020-07-10 [1] RSPM (R 4.0.2)
##  pkgconfig      2.0.3    2019-09-22 [1] RSPM (R 4.0.0)
##  purrr        * 0.3.4    2020-04-17 [1] RSPM (R 4.0.0)
##  R6             2.4.1    2019-11-12 [1] RSPM (R 4.0.0)
##  RBGL           1.66.0   2020-10-27 [1] Bioconductor  
##  rbmn         * 0.9-3    2020-02-13 [1] RSPM (R 4.0.0)
##  Rcpp           1.0.5    2020-07-06 [1] RSPM (R 4.0.2)
##  readr        * 1.4.0    2020-10-05 [1] RSPM (R 4.0.2)
##  readxl         1.3.1    2019-03-13 [1] RSPM (R 4.0.2)
##  reprex         0.3.0    2019-05-16 [1] RSPM (R 4.0.0)
##  Rgraphviz      2.34.0   2020-10-27 [1] Bioconductor  
##  rlang          0.4.8    2020-10-08 [1] RSPM (R 4.0.2)
##  rmarkdown      2.4      2020-09-30 [1] RSPM (R 4.0.2)
##  rmdformats     0.3.7    2020-03-11 [1] RSPM (R 4.0.0)
##  rsconnect      0.8.16   2019-12-13 [1] RSPM (R 4.0.0)
##  rstudioapi     0.11     2020-02-07 [1] RSPM (R 4.0.0)
##  rvest          0.3.6    2020-07-25 [1] RSPM (R 4.0.2)
##  scales       * 1.1.1    2020-05-11 [1] RSPM (R 4.0.0)
##  sessioninfo    1.1.1    2018-11-05 [1] RSPM (R 4.0.0)
##  stringi        1.5.3    2020-09-09 [1] RSPM (R 4.0.2)
##  stringr      * 1.4.0    2019-02-10 [1] RSPM (R 4.0.0)
##  tibble       * 3.0.4    2020-10-12 [1] RSPM (R 4.0.2)
##  tidyr        * 1.1.2    2020-08-27 [1] RSPM (R 4.0.2)
##  tidyselect     1.1.0    2020-05-11 [1] RSPM (R 4.0.0)
##  tidyverse    * 1.3.0    2019-11-21 [1] RSPM (R 4.0.0)
##  utf8           1.1.4    2018-05-24 [1] RSPM (R 4.0.0)
##  vctrs          0.3.4    2020-08-29 [1] RSPM (R 4.0.2)
##  withr          2.3.0    2020-09-22 [1] RSPM (R 4.0.2)
##  xfun           0.18     2020-09-29 [1] RSPM (R 4.0.2)
##  xml2           1.3.2    2020-04-23 [1] RSPM (R 4.0.0)
##  yaml           2.2.1    2020-02-01 [1] RSPM (R 4.0.0)
## 
## [1] /usr/local/lib/R/site-library
## [2] /usr/local/lib/R/library